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ABSTRACT 

We consider the acceleration of very small dust grains including Polycyclic Aromatic Hydrocarbons 
(PAHs) arising from the electrostatic interactions of dust grains that have charge fluctuations in time 
due to charging events. We simulate the charge fluctuations of very small grains due to their sticking 
collisions with electrons and ions in plasma, and the emission of photoelectrons by UV photons using 
Monte Carlo method. We identify the acceleration induced by the charge fluctuations as the dominant 
acceleration mechanism of very small grains in the diff'use interstellar medium (ISM). We show that 
this acceleration mechanism is efficient for environments with low degree of ionization (i.e., large 
Debye length), where the charge fluctuations arc slow but have a large amplitude. We also discuss 
the implications of the present mechanism for grain coagulation and shattering in the diffuse ISM, 
molecular clouds, and protoplanetary disks. 

Subject headings: dust, extinction - ISM: kinematics and dynamics, acceleration - ISM 



1. INTRODUCTION 

Dust is an important constituent of the interstellar 
medium (ISM), molecular clouds, and accretion disks 
(see Whittet 2003; Drainc 2009, 2011). It gets involved in 
many key processes, for instance, it controls heating and 
cooling of the ISM (see Draine 2003; Tielens 2005), re- 
veals magnetic flelds through grain alignment (see Lazar- 
ian 2007 for a review), and interferes with attempts to 
measure the temperature anisotropy and polarization of 
the cosmic microwave background (CMB) radiation (see 
Lazarian & Finkbeiner 2003; Fraisse et al. 2009; Dunkley 
et al. 2009). 

Very small dust grains of size a < 100 A with a notable 
fraction of polycyclic aromatic hydrocarbons (PAHs, 
hereafter very small grains and PAHs are used inter- 
changeably) radiate electric dipole emission, which is an 
important component of CMB foregrounds (Draine & 
Lazarian 1998; Hoang et al. 2010, 2011). PAHs in proto- 
planetary disks (PPDs) around young intermediate mass 
(Hcrbig Ac/Be) stars (Acke & van den Ancker 2004) and 
in outer layers of low mass (T-Tauris) stars (Gcers et al. 
2006; Oliveira et al. 2010) may affect the magnetoro- 
tational instability (MRI) due to their influence on the 
environment ionization (Bai 2011; Perez-Becker & Chi- 
ang 2011). 

Most properties of dust, including light extinction, 
electron photoemission, and chemical activity depend 
not only on grain chemical composition, but also on their 
sizes. In astrophysical environments, grain size is affected 
by grain-grain collisions, which depends on their relative 
velocities. The minimal velocity of grains corresponds 
to their thermal velocity at the ambient gas tempera- 
ture arising from Brownian motions. These velocities 
are usually assumed in the models of dust coagulation 
(Ossenkopf 1993; Dominik & DuUemond 2008). It has 
been long known that large scale hydrodynamic motions 
associated with turbulence can make grains move faster 
(see Draine 1985), but detailed calculations applicable to 
astrophysical environments started to appear only a few 
years ago. 



Since most astrophysical media are magnetized and 
grains are charged, the hydrodynamic treatment of accel- 
eration is not adequate. A proper treatment of the grain 
acceleration through the resonant interactions of charged 
grains with magnetohydrodynamic (MHD) turbulence 
has been developed relatively recently (Lazarian & Yan 
2002; Yan & Lazarian 2003; Yan ct al. 2004; Yan 2009; 
Hoang et al. 2011b). This treatment makes an extensive 
use of the advances in the theory/numerical studies of 
compressible MHD turbulence (see Cho & Lazarian 2003; 
Kowal & Lazarian 2010 and ref. therein) and provides 
the mathematical formalism of the second-order Fermi 
acceleration of charged grains interacting with MHD tur- 
bulence. 

The acceleration due to the resonant interactions of 
MHD turbulence with charged grains decreases with de- 
creasing grain size. Such decrease arises from the fact 
that the Larmor radius of charged grains becomes smaller 
as the grain mass decreases, and correspondingly, grains 
have to interact with smaller and less powerful turbulent 
fluctuations. In addition, compressible fluctuations (i.e., 
fast modes), which were identifled in Yan & Lazarian 
(2003) as the most efficient source of acceleration, get 
suppressed at small scales due to plasma viscous damp- 
ing, while the Alfvcnic mode gets inefficient for acceler- 
ation at the small scales due to anisotropy (see Yan & 
Lazarian 2003 for more discussion). The resonant accel- 
eration for grains with sizes < 10~^ cm becomes rather 
inefficient in most media discussed in Yan & Lazarian 
(2003) and Yan et al. (2004, hereafter YLD04). These 
conclusions were also confirmed by Hoang et al. (2011b) 
who considered the non-linear effects of resonant acceler- 
ation (Yan & Lazarian 2008) and the transit time damp- 
ing acceleration (TTD) of charged grains. 

A new mechanism of astrophysical grain acceleration 
was sketched in Ivlcv et al. (2010), where rough es- 
timates of grain velocities arising from charge fluctua- 
tions in the ISM were provided. The study employed 
Fokker-Planck equation and assumed that (1) charge 
fluctuations are fast, i.e., the relaxation timescale for 
charge fluctuations is much shorter than the gaseous 
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damping time and the Larmor precession period, and (2) 
charge fluctuations have small amplitude (i.e., 6Q <^ Qo) 
such that these fluctuations can be approximated as be- 
ing continuous and modeled by a Gaussian distribution. 
These assumptions are usually not realistic for very small 
grains for which the charge fluctuations are substantial. 
Therefore, the unrealistically high velocities of very small 
grains were obtained. This motivates our present study 
which is intended to determining the actual velocities of 
PAHs in typical astrophysical conditions. 

In this paper, we aim to quantify the acceleration of 
very small dust grains induced by grain charge fluctua- 
tions using numerical simulations. The structure of the 
paper is as follows. We first discuss major processes re- 
sponsible for grain charging and charge fluctuations in 
§2. In §3 we present an algorithm to simulate the charge 
fluctuations using Monte Carlo method. In §4 we present 
a numerical simulation method to investigate grain accel- 
eration, including dynamical equations and Monte-Carlo 
simulations. Velocities of very small grains in the ISM 
from numerical simulations are presented in §5. Discus- 
sions and summary are presented in §6 and 7, respec- 
tively. 

2. GRAIN CHARGING AND CHARGE FLUCTUATIONS 

Charging processes for a dust grain in the ISM consist 
of its sticking collisions with charged particles in plasma 
(Drainc & Sutin 1985) and photo-emission induced by 
UV photons (Drainc & Weingartncr 2001). In the former 
case, the grain acquires charge by capturing electrons and 
ions from the plasma, while, in the latter case, the grain 
looses charge by emitting photoelectrons. After a suffi- 
cient long time, these processes result in statistical equi- 
librium of ionization, and the grain has a mean charge, 
denoted by (Z). The grain charge fiuctuates around the 
equilibrium value (Z). 

To characterize the fiuctuations of grain charge, we 
are interested in the charge distribution function fz{Z), 
which describes the probability of finding the grain with a 
charge Ze.^ Detailed calculations for fz arc presented in 
Draine & Sutin (1985) and Weingartncr & Drainc (2001), 
here we describe them briefly. 

Basically, in a steady state approximation, the distri- 
bution function fz of grain charge is constrained by the 
statistical equilibrium 

fz{Z) [Jp,{Z) + .hon{Z)] = fz{Z + l)Je{Z + 1), (1) 

which means, the number of positively charged particles 
the grain absorbs per second to change the charge state 
from Z to Z + 1 must be equal to the number of electrons 
that the grain absorbs per second to cascade from Z+1 to 
Z. Above Je and Jion are the rate of sticking collisions 
with electron and ion, and Jpe is the rate of emission 
of photoelectron induced by UV photons (see Appendix 
A). To find fz, we solve Equation (1) using an iterative 
method as in Draine & Sutin (1987). 

Let us define a characteristic relaxation time of the 
charge fluctuations, tz, which is equal to the time re- 
quired for the grain charge to relax from Z to the equi- 

^ This distribution function also represents the fraction of grains 
in a charge state because the grain density in the ISM is so low 
that each grain does not interfere with the charging processes of 
the neighboring grains. 



librium state (Drainc & Lazarian 1998b): 



Tz = 



{{Z-{Z)f) ^ al 
J2z fzJtotiZ) J2z fzJtot{Z)' 



(2) 



where Jtot(Z) — Je + Jion + Jpe is the total charging rate 
at the charge state Z. Here we averaged over all possible 
charge states Z to find tz- 



TABLE 1 

Idealized Environments For Interstellar Matter 



Parameters CNM WNM WIM 



riH (cm ^) 

Tgas (K) 

X 

xm 



30 

100 

1 

0.0012 
0.0003 



0.4 

6000 

1 

0.1 

0.0003 



0.1 

8000 
1 

0.99 
0.001 



Here njj is the gas density, Tgas is the gas temperature, xjj is the 
hydrogen ionization fraction, xyi is the metal ionization fraction, x 
is the ratio of radiation energy density to the interstellar radiation 
energy density. 

Figure 1 shows the grain mean charge {the upper 
plot) and charge dispersion dz = <Jzl{Z) {the lower 
plot) as a function of the grain size for the cold neu- 
tral medium (CNM), warm neutral medium (WNM) and 
warm ionized medium (WIM). Physical parameters for 
these phases, including the gas density nn, gas temper- 
ature Tgas, radiation field x, hydrogen ionization frac- 
tion xh and metal ionization fraction xm are given in 
Table 1. The gray area indicates the region where the 
charge fluctuations are important with ct^ > 1. Let 
Ccf be the grain size corresponding to dz = 1- Figure 
1 shows that charge fluctuations can be important for 
a < flcf, where Ocf ~ 10"'^ cm for the CNM and WIM, 
and flcf - 10"^ cm for the WNM. 

We calculate the relaxation time of charge fluctuations 
Tz and the gas damping time Tdrag (Eq. B2) for both 
graphite and silicate grains in various phases of the ISM 
with physical parameters listed in Table 1. Obtained re- 
sults are shown in Figure 2. In the WNM and WIM, 
the charge fluctuations are fast (i.e., tz ^ Tjrag) for the 
entire size distribution. In the CNM, the charge fluctua- 
tions are slow, i.e. tz comparable to Tdrag for a < 10~^ 
cm. The slow charge fluctuations arise from the fact that 
the ionization fraction is rather low in the CNM, which 
makes the charging much less frequent than the gas-grain 
collisions. 

3. MONTE CARLO SIMULATIONS OF CHARGE 
FLUCTUATIONS 

For very small grains (e.g., PAHs) under interest in 
the ISM, the charging is infrequent, so that the charge 
fluctuations may not be adequately characterized by the 
steady state distribution fz ■ In addition, what we are in- 
terested in the present paper is the time-dependent grain 
charge, instead of the steady distribution. Thus, in the 
following, we flrst describe an algorithm to simulate the 
grain charge fluctuations using the Monte Carlo method 
and flnd the time-dependent charge. Then, we compare 
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Fig. 1. — Mean grain charge {Z) (upper) and charge dispersion az 
(lower) of graphite grains as a function of grain size a for various 
ISM phases. Strong charge fluctuations are observed for grains 
smaller than ~ 10~^ cm for the CNM and WIM, and a^i ~ 
10"^ cm for the WNM. 
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Fig. 2. — The gas drag timescale r^jrag and the relaxation time of 
charge fluctuations tz as functions of the grain size a for graphite 
grains in the various phases of the ISM. Charge fluctuations are 
slow with Tz < Ttirag for Q < 5 X 10~* cm in the CNM. 

the resulting charge distribution function with the steady 
state charge distribution obtained using the statistical 
equilibrium approximation. 

3.1. Monte Carlo Simulations 

Monte Carlo simulations of the charge fluctuations for 
dust grains in plasma have been studied in Cui &: Goree 
(1994), but they considered only the charge fluctuations 
arising from collisions with electrons and ions in the 
plasma. Here, we consider a general case of charging 
arising from both sticking collisions and photoemission 
induced by UV photons. 

The underlying idea of MC simulations is, that, we 
identify the moment of the charging events that occur 
randomly,^ and assume that the grain charge changes in- 
stantly at the moment of charging. We keep track of the 
grain charge and the time interval between two successive 
charging events. The algorithm for simulating the charge 

^ A charging event can be the absorption of an electron (ion) or 
emission of a photoelectron. 



fluctuations is presented in Figure 3. At the initial step 
i ~ and time t = 0, the grain is assumed to be neutral 
with charge Ze = 0. We calculate the rate of sticking 
collisions with electrons and ions, Je{Z), Jion{Z), and 
the rate of emission of photoelectrons, Jpe{Z). The to- 
tal charging rate reads Jtot{Z) = Je + J ion + Jpe- The 
waiting time to the next i + 1 charging event is a random 
variable drawn from a Gamma distribution function with 
the mean value At^ = 1/ Jtot(,Z). 

In fact, the probability to have a charging event occur- 
ring in [t, i + dt] is 



dP = Ar^^ exp {-t/An) dt. 

The probability to find the next charging event i - 
a time interval Ati becomes 



P = 1 



cxp 



At 

"a^ 



(3) 
- 1 after 

(4) 



To find Ati, we generate a random number R in the range 
[0.1], and assign P = R. From Equation (4) we obtain 



AL; = 



ln(l - R) 



(5) 



When the waiting time to the next charging event Ati 
is known, we have to determine what particle among 
e, ion, and photoelectron (pe) the grain captures or 
emits at the next charging event. We note that the prob- 
ability of capturing or emitting u particle is P^ ~ Ju 
with oj = e, ion and pe. The probability that the next 
charging event corresponds to the absorption/emission of 
the X particle among e, ion and pe is given by Px/Ptot 
where Ptot = Pe + Pion + Ppe ■ The next charging event 
corresponds to the absorption/emission of the x particle 
when Px is maximum. 

We begin with generating a random number i?i in 
the range [0, 1] from a uniform distribution. If Ri < 
Praaxl Ptoti thcn wc know that a charging event corre- 
sponding to Pmax just OCCUrrcd. Now, if Pmax = Pe, 

then one electron is captured by the grain in this charg- 
ing event, and the grain charge Z is advanced by —1. If 
^niax 7^ Pe , then the grain captures either an ion or emits 
a photoelectron in this charging event. For this case, Z 
is advanced by -1-1. 

If Ri > Pmax/ftot, then we know that the next charg- 
ing event with Pmax has not occurred, and the charg- 
ing event occurred with one of two remaining parti- 
cles. To determine what particle the grain actually ab- 
sorbs/emits, we continue by generating a random num- 
ber i?2, and compare it with P of these particles. The 
grain charge Z is updated, and the time is advanced as 
ii+i = ti + Ati. We tabulate the charging time ti+i 
and grain charge Z{t). This process is iterated until the 
elapsed time t^+i equal to a few gaseous damping times 

'^drag ■ 

Figure 4 shows the time-dependent grain charge ob- 
tained usin g MC simulations for a 10, 20 and 100 A gram 
in the WIM. The discrete nature can be easily seen for 
very small grain a = 10 A in which Z varies slowly. 

3.2. Charge distribution function 

To characterize the fluctuations of grain charge from 
MC simulations, we find the steady state charge distri- 
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Fig. 3. — Diagram for MC simulations of charge fluctuations. 
Initially, grain is neutral. The charging rates Je, J ion a^nd Jpe 
are used to calculate the mean time between two charging events 
Ar and the waiting time to the next charging event Af. The 
simulations are iterated until the elapsed time tij^\ exceeds the 
drag time tdrag- 



bution function jz using the tabulated data Z(t). Fig- 
ure 5 compares obtained from MC simulations with 
jz obtained using the statistical equilibrium approxima- 
tion (see Eq. 1) for different grain sizes in the ISM. 
For a = 100 A grains, the results from two different ap- 
proaches are similar, while there exists some difference 
in ]z for smaller grains (e.g., a = 10 and 20 A). Such 
a difference stems from the fact that the charge fluctu- 
ations occur at much lower rate for very small grains, 
so that the statistical equilibrium approximation is not 
adequate. 

4. NUMERICAL SIMULATIONS OF GRAIN 
ACCELERATION 

4.1. Basic Equations 

Let consider an ensemble of grains in a plasma with 
density uh, ionization fraction xh, and temperature Tgas. 
The velocity of a grain i of mass changes due to the 
gas drag force as well as the Coulomb force induced by 
charged grains in its vicinity. The conventional equation 
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Fig. 4. — Grain charge as a function of time obtained from MC 
simulations for three values of grain size a = lOA (red line), 20 A 
(blue line), and 100 A (black line) in the WIM. Charge fluctuations 
are rather slow for the smallest grains (red and blue lines) , but fast 
for larger grain (black line). 
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Fig. 5. — Charge distribution function fz obtained from MC 
simulations (red line) compared to fz from the statistical equilib- 
rium approximation (black line) for three values of grain size in the 
WIM. fz from two approaches is almost similar for the 100 A grain 
when the charge fluctuations are fast, but there is some difference 
for two smaller sizes for which the charge fluctuations arc slow. 



of motion reads 



Tdr, 



F 

m,. 



(6) 



where Tjiag is the damping time due to gas drag (see 
Eq.B2), Ri is a random force due to Brownian motions, 
which is described by 



{Ri) - Gi, 



(7) 



with Gi = 2rgasfcB?gas being the translational excita- 
tion coefficient due to dust-gas collisions. Here is the 
Coulomb force acting on i-grain induced by all surround- 
ing grains j ^ i- Assuming the Yukawa potential with 
the screened length A, F^ is given by 



Qi{t)Qj{t)e^p {-nj/X) 



5 



exp(-r,;j/A) 



(8) 



where Qi(t) and Qj{t) are the charge of i and j grams 
at the time t, and r.y is the separation between them. 

The position of the i grain is determined by the equa- 
tion 

Solving Equation (6) is rather challenging, because it 
is involved with the stochastic term Ri and a time- 
dependent force F, which depends on the instant grain 
charge Q{t). 

4.2. Numerical Integration 

To solve Equation (6) combined with Equation (8) for 
grain velocities, we first need to know the fluctuating 
charge Q{t) at the time t. Detailed treatment of Q{t) 
using data from MC simulations will be addressed in 
the next section. Here, we present a general numerical 
approach to solve Equation (6) for an arbitrary time- 
dependent force F(r, i). 

Let us denote /i„ be the time interval at the step n in 
which the charge of all grains are constant. The velocity 
of a grain at the step n -I- 1 is given by 



f n+l — r„ + V„/l„, 



'n+l 



tn+1 



Ri 



2m, 



1 ( 



where F^ = Fi(r„,i„) and F[ = F,(f„+i,t„). 
In simulations we adopt dimensionless units 

t , V 

, V , 

Tdrag Vt 



t' 



(10) 
(11) 

(12) 



(13) 



1/2 

where vt = (fceTgas/'Tii) is the grain thermal velocity. 
The position is normalized over the screened length scale 
A: 



(14) 



Using the dimensionless units. Equations (10)-(12) can 
be rewritten as 



•■n+l ■ 



'n+l 



■■n+l 



R' 



'^drag 



2mi dtA^ 



h' 



hi 



(15) 
(16) 

(17) 



When the velocities of grains is known, the kinetic tem- 
perature is evaluated by the ensemble average: 



1 ^ z 

1 X - rudV 



N ^ 



(18) 



In the above equations, we note that the grain charge 
Q changes in time. Therefore, to obtain grain velocity at 
each timestep, we need to treat Q{t) properly. Below, we 
present an algorithm to find Q(t) at each timestep using 
data from MC simulations. 



4.3. Algorithm 
4.3.1. Tiny grains a < 10 A 

Very small grains with a < 10 A have infrequent charg- 
ing (see Fig. 4), thus most of the time they spend in the 
charge states Z = 0,±1, and ±2. 

Since the grain charge Q is a function of time, the time 
step hn in Equation (15)-(17) must be chosen such that 
the charge of grains in the simulation box remains con- 
stant during this timestep. The easy way is to choose the 
timestep coincident with the charging time of all grains. 

We consider Np grains in the simulation box. For each 
grain, wc generate Nz random charging events, and tab- 
ulate the moment of charging tz and the time interval 
between two charging events Atz using MC simulations 
(here tz = ti in Section 3). As a result, we have total 
-^chrg = Np X Nz charging events in the entire simu- 
lation box. When the charging moment tz is known, 
we can sort it in ascending order, and obtain f^^^j-g with 
n ~ 1 — A^chrg for the sequence of A^chig charging events in 
the simulations. The timestep is then chosen as the time 
interval between two successive charging events in the 
sequence, i.e., /i„ = t^^^ - 1"^^^^. When sorting ichrg, we 
also mark the grain that experiences the charging event. 
Therefore, at each time step, only this marked grain has 
charge that changes instantly by ±1, while the charge 
of other grains remains the same. When the integration 
time t„ is equal to the charging time tz,i, then the charge 
of i grain changes. 

If the timestep hn results in a large increment of d„, 
then we need to adjust /i„ so that the increment of grain 
velocity is not larger than an e fraction of Vn ■ We choose 
e = 0.1 for our simulations. When the time interval 
between charging events is very large, we need to limit 
h'^^^ = 10~^ to keep the numerical integration stable. 

4.3.2. Very small grains a > lOA 

For grains larger than 10 A, the charging occurs more 
frequently, so that the time interval between two succes- 
sive charging events dtchrg may be rather small (about 
lO^^Tdrag)- As a result, we need to integrate Equation 
(6) over a huge number of timesteps to achieve the termi- 
nal velocity after Tdrag, which is rather impractical. To 
overcome this obstacle, we choose a constant timestep 
hn, and employ the MC method to generate the instant 
charge at each timestep. 

The basic idea is that, the instant charge Z is ran- 
domly distributed, and can be drawn from the steady 
state distribution fz- We first calculate the accumula- 
tive probability of finding the grain with charge < Ze 
as 



PiZ)= fz{Z.), 



(19) 



where ^^min is the minimum grain charge. The grain 
charge Z at each timestep is random, and obtained by 
solving equation P{Z) = R where i? is a random variable 
drawn from a uniform distribution. 

4.4. Boundary and initial conditions 

We adopt the periodic boundary conditions for our 
simulations. The periodic boundary assures that when a 
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dust grain leaves the box at one boundary, another grain 
is added to the box at its opposite boundary. The peri- 
odic boundary is useful for simulating large scale systems 
consisting of a huge number of particles. 

We start simulations with Np particles randomly dis- 
tributed in a square box of L^. Their initial velocities 
are drawn from a Gaussian distribution with a kinetic 
temperature Tq. 

4.5. Code benchmarking: Brownian motion 

We test the code by simulating the Brownian motions 
of dust grains in a plasma of temperature Tgas- We con- 
sider an ensemble of TV neutral grains of equal size in 
the conditions of the CNM and WIM. Simulations are 
started with an initial temperature Tq < Tgas- We found 
that the grain kinetic energy is driven to thermal equilib- 
rium with the plasma after about the gaseous damping 

time Tdrag- 

5. GRAIN ACCELERATION DUE TO CHARGE 
FLUCTUATIONS IN THE ISM 

5.1. Model parameters 

Below we are going to find the actual velocities of very 
small grains induced by charge fluctuations in various 
phases of the ISM, which have physical parameters listed 
in Table 1 . Wc consider the translational damping due to 
grain-neutral and grain-ion collisions only (see Appendix 
B), and disregard subdominant damping processes (e.g., 
dipole-dipole, grain-dipole interactions). We also disre- 
gard the sticking grain-grain collisions due to the fact 
that the mean distance between grains in the ISM is 
large, which is comparable to the Debye length (see Ivlev 
et al. 2010). 

5.2. Grain size distribution 

First, let us assume for the sake of simplicity that all 
dust grains have the same equivalent size Oeq with density 
'T'd(aoq)- Here, the density nd(aeq) is constrained by the 
condition that the mass density of acq grains is equal to 
the total mass density of dust integrated over the entire 
size distribution from amin to a,„ax: 

47r o , , /■°""'- 47r , 
A/d = Y<qnd(aeq) = j Y d^ ' ^ ^ 

where flmin = 3.56 A and Omax = 2.5 x 10^^ cm are 
assumed. Assuming the MRN distribution dn/da = 
JT-H^MRNfl"^'^ with Amrn = 10~^'^ -^^cm~^-^ for diffuse 
clouds (see Mathis et al. 1977), we obtain 

nd(acq)«10-V(T7^)~'. (21) 
V 10 'cm/ 

We perform MC simulations for different equivalent 
sizes Coq from Omin to Cmax- For each value of acq, the 
density is inferred from Equation (21). 

Next, we account for the electrostatic interactions be- 
tween grains of different sizes. We divide the grain size 
distribution from to fln^ax into iVbin with a constant 
log binsize, Alna = ln(a,nax/amin)/(A^f)m)- The grain 
size in the i bin is given by 

Inai = Inomin + iAlna, (22) 



and the fraction of grains in i bin to the total dust density 
Ud reads 

ri£ ^ {dn.i/da)Aai _ a'^'^aj , . 

nd rid ^^^^ a^^ Ok 

where we have adopted the MRN distribution for sim- 
plicity. 

Assuming the total number of particles in the simula- 
tion box is Np^ the box size L is constrained by 

A^bin-l 

2 n,L^ = UdL" = Np. (24) 

We run simulations for Np grains of various size a, and 
find the average velocity for grains having the same size. 
For each grain size a, the dominant contribution to the 
acceleration arises from its interaction with other grains 
of size a' > a. Therefore, the Coulomb force in Equation 
(8) corresponds to the summation over all grains larger 
than the grain of mass mi under interest. 

5.3. Grain Velocities 

We present the results of MC simulations for first the 
simplest case in which all grains have the same equivalent 
size. For each grain size, we follow the time evolution of 
grain velocity v and calculate the averaged velocity v at 
each time step. We identify the time interval in which v 
varies slowly and then saturates. The terminal velocity 
is obtained by averaging v over the time interval ranging 
from Tdrag to the total integration time T . Here, we adopt 
the total integration time T = 2Tdrag- 

Figure 6 presents the temporal evolution of grain veloc- 
ity normalized over its thermal velocity in the CNM and 
WIM for three grain sizes. As expected, the grain veloc- 
ity increases exponentially as a function of time from the 
initial value, then slows down gradually and saturates af- 
ter about a gas damping time. The terminal velocity is 
obviously larger for smaller grains. In addition, we can 
see that the acceleration by charge fluctuations is more 
efficient in the CNM than the WIM, which can accelerate 
grains up to about four times its thermal velocity. The 
reason for that is as follows. First, the CNM has stronger 
charge fluctuations, which corresponds to larger charge 
dispersion <tzI {Z) (sec Fig. 1), resulting in stronger fluc- 
tuating electric flclds. Second, since the CNM has lower 
degree of ionization, which corresponds to larger Debye 
length, the shielding effect of plasma on electrostatic in- 
teractions between charged grains is much less than in 
the WIM. 

To see the dependence of the terminal velocity with 
grain size a clearly, in Figure 7, we show results for three 
phases of the ISM. We truncate at the size a = 10~^ cm 
because the charge fluctuations are negligible for large 
grains. As shown is the rapid increase of velocity with 
decreasing grain size. 

To study the effect of grain size distribution on grain 
acceleration, we perform simulations for 16 grain size- 
bins from a = 3.56 A to 100 A and take into account 
the electrostatic interactions of grains of different sizes. 
Obtained velocities for the CNM, WNM, and WIM are 
shown in Figure 8. As shown, the velocities of the small- 
est grains increase while the velocities of larger grains 
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Fig. 6. — Evolution of grain velocity normalized to its ther- 
mal velocity as a function of time from MC simulations for three 
grain sizes in the CNM and WIM. Graphite grains are considered. 
Charge fluctuations accelerate very small grains in the CNM to 
suprathermal velocities after about the gas drag timescale r^rag- 
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Fig. 7. — Terminal grain velocities as a function of grain size for 
the CNM, WNM and WIM due to charge fluctuations arc shown 
by solid, dotted and dashed lines, respectively. Thin lines show the 
thermal velocities of grains in the corresponding ISM phases. 

decrease slightly, compared to the velocities obtained for 
the grains of equal size. That is understandable because, 
when the grain size distribution is taken into account, 
the smallest dust grains (PAHs) have chance to interact 
with larger grains. Largest grains, on the contrary, in- 
teract mostly with smaller grains so that their velocities 
decrease. 

5.4. Comparison with Fokker-Planck Equation 
Approach 

The Fokker-Planck (FP) equation was used in Ivlev et 
al. (2010) to study the grain acceleration by charge fluc- 
tuations. The key assumptions in the approach based 
on the FP equation are as follows: (1) grain charge fluc- 
tuations are fast and (2) the charge dispersion is small 
compared to the equilibrium charge. These assumptions 
are, in general, valid for large grains, but they fail for 
very small grains (see Fig. 1). To see the consequence of 
such unrealistic assumptions on predicting PAH veloci- 
ties, we conservatively apply the FP approach for very 
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Fig. 8. — Similar to Fig 7, but the effects of grain size distribution 
are taken into account. Grain velocities predicted by the Fokker- 
Planck equation are much larger than the actual values computed 
by MC simulations. 

small grains and compute grain velocities using Equation 
(CIS) in which the actual values of az and (Z) from Fig- 
ure 1 are adopted. 

Figure 8 compares grain velocities obtained from MC 
simulations with those predicted by the FP approach for 
the ISM. We can see that the FP approach predicts much 
larger grain velocities than the actual values from MC 
simulations. The difference in grain velocities from the 
two approaches decreases when a increases and becomes 
negligible for sufficient large grains. 

The reason for that the FP approach predicts unreal- 
istically high velocities for PAHs is the following. When 
adopting the FP approach to estimate grain velocities, we 
are implicitly assuming an unrealistic situation in which 
PAHs have fast fluctuating charges with large amplitude. 
It corresponds to the situation in which each PAH is 
subject to unusually strong and fast fluctuating electric 
fields. As a result, PAHs gain a substantial amount of 
energy after a short interval of time and are accelerated 
to extremely high velocities. However, for very small 
grains (e.g., PAHs), we showed in Figures 1 and 2 that 
the grain charge indeed has large dispersion, but it fluc- 
tuates rather slowly, which makes the acceleration much 
less efficient than the prediction by the FP equation. 

For highly ionized media, the fast charge fluctuations 
can be applied, but we found that the results are still 
much lower than the FP theoretical prediction. Such a 
difference can arise mainly from the fact that the plasma 
screening effect is disregarded in the estimate of Ivlev 
et al. (2010), while the highly ionized media have short 
Debye lengths, which reduce the electrostatic dust-dust 
interactions. 

6. DISCUSSION 

6.1. Grain Acceleration from Various Mechanisms 

Grain acceleration by hydrodrag (see Draine 1985) was 
discussed as the dominant mechanism driving grain mo- 
tions in the ISM as well as in protoplanetary disks. Re- 
cently, thanks to signiflcant progresses in understanding 
of MHD turbulence, Lazarian & Yan (2003), YLD04, and 
Yan (2009) identified the gyro-resonant interactions of 
grains with fast modes in MHD turbulence as the dom- 
inant mechanisms to accelerate large grains (a > 10~^ 
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cm) to super- Alfv'enic velocities. For very small grains 
(e.g., PAHs and nanoparticles) , the gyro-resonance accel- 
eration is inefficient for most phases of the ISM, because 
PAHs and nanoparticles with smaller gyro radii only res- 
onantly interact with small scale and weak fluctuations. 
Hoang et al. (2011b) revisited the treatment of gyroreso- 
nance acceleration by MHD turbulence by accounting for 
the fluctuations of the grain guiding center with respect 
to the mean magnetic flcld. They improved estimates 
of the previous authors, and proposed a new way of ac- 
celeration through transit time damping (TTD) of fast 
modes, which is efficient for super- Alfvenic grains (see 
solid lines in Fig. 9). 

6.2. Acceleration induced by Charge Fluctuations 

In the present paper, we have numerically investigated 
the grain acceleration mechanism induced by charge fluc- 
tuations for very small grains. The combination of MC 
simulations for charge fluctuations with direct simula- 
tions of electrostatic interactions between grains with 
fluctuating charge allows us to follow the evolution of 
grain velocities in time and estimate their terminal ve- 
locities. We showed that the acceleration induced by 
charge fluctuations is indeed an important mechanism to 
accelerate PAHs. Resulting grain velocities can be sev- 
eral times higher than their thermal velocities, which is 
obviously important. 
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Fig. 9. — Grain velocity as a function of grain size a due to charge 
fluctuations (shaded area) and due to resonance acceleration by 
MHD turbulence (Hoang et al. 2011b) for the CNM, WNM, and 
WIM. MHD turbulence is the efficient mechanism for large grains, 
and grain charge fluctuations work efficiently at the small end of 
grain size. Dotted and solid lines correspond to the acceleration 
from gyroresonance and gyroresonance plus transit time damping. 

We also found that grains velocities obtained from MC 
simulations are between two and three orders of mag- 
nitude lower than the prediction by the FP equation. 
Such overestimate of the FP approach arises because the 
charge fluctuations of very small grains in the ISM are 
slow, in contrast to the assumption of fast charge fluctua- 
tions assumed in the FP treatment, which makes electric 
fields fluctuating rather slowly. In addition, the FP treat- 
ment disregarded the screening effect of plasma. In the 
ISM, the mean separation between two grains is compa- 
rable to the Debye length. Thus, the Coulomb interac- 
tion is suppressed considerably, and the grain accelera- 
tion is decreased accordingly. 



Figure 9 compares grain velocities induced by charge 
fluctuations (shaded area) with those due the resonance 
interactions in MHD turbulence. As shown, MHD turbu- 
lence is efficient for large grains, while charge fiuctuations 
are important for very small grains. 

6.3. Grain coagulation and shattering 

The grain coagulation is a process in which small dust 
particles collide with each other to form aggregates of 
submicron size (Spitzer 1978). The grain coagulation 
can occur in dense clouds (e.g., Ossenkopf 1993; Ormel 
et al. 2009) as well as in protoplanetary disks (see van 
Boekel et al. 2003). 

Grain coagulation and shattering result from grain- 
grain collisions, which depends on grain relative veloc- 
ity, denoted by Vdd- The threshold velocity for the grain 
shattering is a function of the grain size: 

«shat=2.7 — ^ kms-i, (25) 

V 10~' cm/ 

(Chokshi et al. 1993). If Vdd < Wshat, the grains col- 
lide and stick together. When Vdd > Wshat, the collisions 
with high velocity produce shock waves inside the grain, 
and shatter them in smaller fragments. For Vdd — ^ 20 
kms~^, the evaporation of dust grains occurs and grains 
are destroyed. 

We note that due to resonant acceleration by MHD 
turbulence, large grains a > 10~^ cm move frequently 
with velocity larger than Wghat (see Fig. 9). Thus, the 
shattering of large grains can be efficient to replenish very 
small grains into the ISM. The effect of grain acceleration 
by MHD turbulence on grain coagulation and shattering 
was studied in Hirashita & Yan (2009). In some phases, 
they found that grain acceleration can modify grain size 
distribution. At the same time, very small grains move 
faster than thermal speed due to the acceleration dis- 
cussed in the present paper and enhance grain coagula- 
tion. Therefore, future models of grain evolution should 
take into account the acceleration of very small grains. 

6.4. Effects of PAH acceleration in protoplanetary disks 

Grain coagulation is the first step to planetesimal for- 
mation in protoplanetary disks (PPDs). It is also be- 
lieved that the coagulation from very small grains to mi- 
cron size occurs rapidly (see e.g., van Boekel et al. 2003), 
so that PAHs and nanoparticles are rapidly depleted in 
PPDs. However, recent observations by Infrared Space 
Observatory (ISO) and 5piteer satellite reveal PAH emis- 
sion features, which indicates the existence of PAHs in 
young intermediate mass (Herbig Ae/Be) stars (Acke & 
van den Ancker 2004) and in the surface layers exposed to 
UV radiation in young low mass (T-Tauri) stars (Geers 
et al. 2006; Oliveira et al. 2010; Berne et al. 2009). 

PPDs have highly stratified structure, that means, 
there exists some gradient of grain size from the mid- 
plane to the surface. When PAH emission is seen, it 
is probably only from PAHs in the surface layers that 
"see" direct stellar radiation. The abundance of PAHs 
and other nanoparticles below the surface is rather un- 
certain. 

(a) Effects of acceleration on dust coagulation and 
planet formation 

Current models of grain growth in PPDs disregarded 
grain charges and their fluctuations. Recently, Okuzumi 
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(2009) took into account the effects of grain charge 
and showed that Coulomb potential barrier can suppress 
grain coagulation. Okuzumi (2009) appealed to turbu- 
lence as an energy source to overcome the potential bar- 
rier. Below, we show that charge fluctuations can induce 
very small grains to overcome the Coulomb barrier as- 
suming that PPDs have a sufficiently high degree of ion- 
ization for which the charge fluctuations are important. 

Assuming that the charge fluctuations are Gaussian, 
then, the mean charge and charge dispersion are respec- 
tively given by 



-0.1 



10- 



100 K 



(26) 



where z = 2.5 is adopted (see MorfiU & Ivlev 2009). 

The dust mass density in PPDs is related to the gas 
density as 



Pd — fdgPg, 



(27) 



where fdg = 0.014 for the PPD (see Tanaka et al. 2005). 
The dust density of the grain size a is then 



na - 10 



V 10-6 cm/ V; 



(28) 



Using above equations for Equation (C16), we obtain 

2 



_ n oi -|nl7 2 ( ^gas Y ( a_ ^ ^ 

-3X 2 



X Xi 



0.1 cm 



cm 



1010 cni- 



(29) 



where e < 1 is introduced to account for the overestimate 
of grain velocities obtained using the FP approach, Ui is 
the ion density, Xi = ni/rigns is the ionization fraction 
of gas, and PPDs are assumed to be isothermal plasma 



with Tgas = Ti. 



Now, let us consider whether the charge fluctuations 
can help grain grow against the Coulomb potential bar- 
rier. The ratio of the grain kinetic energy arising from 
charge fluctuations to electrostatic potential at a distance 
equal to 2a with a being the grain size is 



From Equations (26), (29), and (30) we obtain 



(30) 



i^Coul 



^ 2.3 X lO^e^ f ^ ) 

ui V 10-6 cm/ 

lOio cm-3 



100 K 

10-4 



Xi 



(31) 



where the ionization fraction Xi was assumed to 
be sufficiently high so that the charge fluctuations are 
important. 

Grain growth occurs only if £'kin/-£'Coui 1^ 1- Thus, 
from Equation (31) we can derive the range of grain size 
with the grain growth. Using parameters for PPDs con- 
ditions, Tgas = 300 K, rigas = lOi" cm-3 and x, = IQ-^ 
and taking e — 10" i, we found that very small grains 
with size a < 6 x 10-^ cm. Therefore, the growth of 



larger grains can be induced by other processes (e.g., 
turbulence as suggested in Okuzumi 2009). 

The presence of PAHs and nanoparticles in PPDs can 
arise from the fragmentation and shattering of micron 
graphite and silicate grains (see e.g., DuUemond & Do- 
minik 2008). The level of MHD turbulence in PPDs is 
very uncertain, but it is believed that the layers near the 
surface are highly turbulent. As a result, large grains can 
be accelerated by the resonant interactions of fast modes 
to super-Alfvenic speed (Yan et al. 2004; Hoang et al. 
2011b). Thus, the grain shattering may be efficient that 
return very small grains to PPDs. 

(b) Effects of PAH acceleration on magnetorotational 
instability 

PAHs play an important role in the ionization equi- 
librium in PPDs. Bai (2011) found that PAHs can fa- 
vor the accretion due to magneto-rotational instability 
(MRI) by decreasing ambipolar diffusions in outer layers 
of PPDs. In fact, if charged PAHs are more abundant 
than electron (npAH > n-e), then the ambipolar diffusion 
is less important. However, since the electrical conduc- 
tivity by charged PAHs is more efficient than electron, 
the MRI becomes more efficient. Suprathermal velocities 
of PAHs due to charge fluctuations in PPDs can enhance 
the accretion rate of electron and ion onto grains, which 
certainly affect the ambipolar diffusion^ and MRI. Thus, 
understanding grain motions in PPDs is important for 
understanding the dynamics of PPDs and planetesimal 
formation. 

7. SUMMARY 

The present paper applies the Monte Carlo simulations 
to quantify the efficiency of the new acceleration mecha- 
nism for very small grains induced by charge fluctuations. 
Our principal results can be summarized as follows. 

1. Grain charge fluctuations due to both the sticking 
collisions with electrons and ions in the plasma and the 
emission of photoelectrons by UV photons are simulated 
using the Monte Carlo method. The charge distribu- 
tion obtained by MC simulations is similar to the steady 
charge distribution for large grains, but they are different 
for very small grains. 

2. Grain acceleration arising from electrostatic in- 
teractions between dust grains is investigated using the 
Monte Carlo simulations. We found that the acceler- 
ation induced by charge fluctuations is efficient for very 
small grains, which can increase grain velocities to several 
times their thermal velocities. This mechanism is more 
efficient for the CNM with lower ionization fraction (i.e., 
larger Debye length) and larger charge dispersion than 
the WNM and WIM. 

3. The increase of relative velocities of PAHs and 
nanoparticles due to charge fluctuations that we reported 
has certain effects on grain coagulation and shattering 
and should be accounted for in the models of grain evo- 
lution and planetesimal formation. 

AL and TH acknowledge the support of the NSF- 
funded Center for Magnetic Self- Organization (CMSO) 

^ The ambipolar diffusion is usually invoked to describe the re- 
moval of magnetic field during star formation. Another powerful 
mechanism is reconnection diffusion based on the ability of tur- 
bulent magnetic field to reconnect and redistribute the entrained 
matter (Lazarian 2005; see also Lazarian 2011). 
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APPENDIX 
A. GRAIN CHARGING 

In plasma, grains can acquire charge due to sticking collisions with electron and ions (see Draine & Sutin 1987, 
DS87), and emissions of photoelectrons due to UV photons. 

Al. Collisional charging 

First, let us briefly describe the collisional charging. Consider an beam of charge particles p with velocity v colliding 
with a grain of charge Q ~ Ze, the number of particles sticking to the grain per second is defined as the collisional 
charging rate: 

Jp,v = SpUpVacoii = SpUpVirbl^, (Al) 

where Sp is the sticking coefficient, Up is the number density of particles, and 6cri is the critical impact factor such 
that the closest distance between the particle and the charged grain is r = a. If ZpZ < 0, then we easily get 

^cii ^ ("'^ ~ "rnv'^^a ) values of velocity v. For ZpZ > 0, we have fecri ^ (-'^ ^ wf» ^ a ) ^'^^ velocity > mw^j./2 = 

ZpZe'^/a, i.e., kinetic energy exceeds the repulsive Coulomb energy. 
Assuming that particles have a Maxwellian distribution, we can obtain the charging rate 

Jp{Z)^ J SpnpVacoiifMw{v)4T^v'^dv, (A2) 

where frnwiv) = (fceTgas/STr) exp [—mv'^/2k^Tgas^ is the Maxwellian distribution, and the integral limit is from 
— oo for Zp < and from Wcii ~ oo for ZpZ > 0. 
As a result, we obtain the charging rate from electron and ion sticking collisions: 



1/2 



Jinn = S,:71.; TTfl 



,2 



1 



/ Ze 



for Z < 0, and 



\akTt 
/ Ze 



(A3) 



fSknTeY^^ ^ ( Ze^ . 
Je = Seng ( -33^ ) 7ra exp ( ) , (A4) 

(A5) 



Jion = s^Ui ( ) 7ra exp ( ^ ) , (A6) 



SksTe 



1/2 



2 



Ze^ 
akTp 



(A7) 
(AS) 

for Z > 0. Here Si is ion sticking coefficient, which is taken to be unity (see DS87), Ti and Tg are ion temperature 
and electron temperature. For isothermal plasma, Ti — Tp — Tgas. Above we ignore the minor contribution from the 
interaction of charge with the image potential on the grain. Detailed descriptions can be found in given in DS87. 

A2. Photoemission 

The number of electrons leaving the grain due to photoemission by UV photons per second is 

/u c 
YeQ^ksj^d,^, (A9) 

where = u^/ {hu) is the density of photon with energy E = hv > 13.66^, and Yp is the photoemission yield, and 
Qahs is the absorption efficiency (see detail in Weingartner & Draine 2001). For the ISM, we calculate photoemission 
rate Jpp using the interstellar radiation from Mathis et al. (1983). We adopt the tabulated data Qabs for graphite 
grains from Draine & Li (2001). 

In the ISM, the photoemission is dominant, but in dense regions with high density (e.g., protoplanetary disks) the 
collisional charging is dominant in mid-plane regions where UV photons from the new born stars are shielded. Close to 
the young stars or higher from the mid-plane, the photoemission can be dominant due to abundance of UV photons. 
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B. GAS DRAG 

Interactions of dust grains with the ambient gas represent the primary mechanism of dissipating streaming motions 
of grains. The damping rate of translational motion arising from the interaction with neutral gas is essentially the 
inverse time for collisions with the mass of the gas equal that of a grain (Purcell 1969), 

where m„, n„, and T„ arc the mass, volume density, and temperature of neutral atoms, is the mass density of dust 
grains, and a is the grain size. 

When the ionization is sufficiently high, the interaction of charged grains with ions in plasma becomes important 
(Draine & Salpeter 1979). The ion-grain cross section due to long-range Coulomb force is larger than the atom-grain 
cross section. As a result, the rate of translational motion damping gets modified (Draine & Salpeter 1979).^ 

For subsonic motions total damping rate due to neutral and ions is 



r. 



.-1 



drag 



Idg = 7" + 7jon = Q!7„, (B2) 



where the rcnormalizing factor reads 



z 



, 3 {Mf/^ 

2V^|Z|e3(xnH)i/2 ' ^^"^^ 

Here Xi is the abundance of ion i (relative to hydrogen) with mass and temperature T^, x = Xi, and f{Z) is 
the grain charge distribution function. When the grain velocity Vd relative to the gas becomes supersonic the dust 
interactions with the plasma is diminished, and the damping rate in this case is renormalized due to the gas-dynamic 
correction (Baines et al. 1965; Purcell 1969; Draine & Salpeter 1979) 



where Cg = ^Jk^Tn/rrin is the sound speed. 

C. THEORETICAL TREATMENT FOR ACCELERATION DUE TO CHARGE FLUCTUATIONS 

CI. Charge evolution for large grains 

In Section 2 and 3, we discussed the charging and charge fluctuations for very small grains and present Monte 
Carlo simulations to model the charge fluctuations. Here, we consider the case of large grains for which the charge 
fluctuations can be modeled using Langevin equations. 

Assume that some small deviation from the mean grain charge is introduced, then the charging rate changes accord- 
ingly to bring the charge to an equilibrium state. Denote Tch be the relaxation time for which the charge fiuctuations 
return to the equilibrium state, the charging frequency becomes 



^ Q=Q„ 



where Jk with fc=ion, e and pe, is the charging rate from the k process. 

For the collisional charging, using Equations (A3) and (A4) for (CI), we obtain Vch and the charge dispersion ctq 
for a grain of size a 

d{Jion - Je) 1 + z a 
Vc\Y^ -J7\ = — 1 — (C-^) 



dQ 2-K Ad 

l + z 
z{2 + z) 



- ^^leQol, (C3) 



where z = Qoe/ (akTi) and 



This drag consists of both close collisions and Coulomb distant r 1 4 ■ 1, 11 j-i, j- t ■ j j i iv ■ 

I , . o- ,1 , 01 electrons IS much smaller than that or 10ns, drag due to collisions 

interactions ot dust with 10ns and electrons, since the momentum ii, 1 -t j- j j > =. 

with electrons are disregarded. 
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is the Debye length due to ion with the plasma frequencies 



pi 1 ^vd 



(C5) 



C2. Continuous Approximation 

The charge fluctuations of large grains can be described by a Gaussian and Markovian random process. The time 
evolution of charge can be modeled by a stochastic differential equation: 

^^-'^chQ + Tq, (C6) 

where 

(r(i)r(t + r)) =2i.ehfTQ(5(r). (C7) 

where ctq is the standard deviation of charge from equilibrium or charge dispersion (see a review by Morfill & Ivlev 
2009). The auto-correlation function of charge is a decaying law 

{Q{t)Q{t + r)) = Qge.Tp(-i.ehT). (C8) 



C3. Evolution of grain kinetic temperature 

In order to understand the effect of the energy variation in binary dust-dust collisions on the mean kinetic energy of 
the whole dust ensemble, one should employ the kinetic approach. The kinetics is described in terms of the velocity 
distribution function fd{p,t). There are two principal contributions to the dust kinetics ~ one is due to the mutual 
dust collisions and another because of dust interactions with the ambient gas. 

Three known interaction processes between dust-gas are dust-neutral, dust- ion collisions, and plasma drag, which 
arises from interaction of passing ions with grain dipole moment. We assume that these dust-gas interactions lead to 
the local thermal equilibrium of gas and dust, characterized by a temperature Tg, at different rates 7n,7i and (see 
Appendix B). 

The resulting kinetic equation has the following form: 

~ Stddfd + Stdgfd, (C9) 

where Stdd and St^g denote the collision operators (integrals) describing the dust-dust and dust-gas interactions, 
respectively. Although the analysis of Eq. (C9) can be performed for arbitrary fd{p), for the sake of convenience we 
assume that dust particles have a Maxwellian velocity distribution /m(p)- Then the equation for evolution of mean 
kinetic energy (kinetic temperature) kTd{t) = J {p'^ /md)fMdp is obtained by taking the second moment of Eq. (C9), 



Td= \ I — {Stddfu + StdgfM) dp. (CIO) 
3 J md 

Ivlev et al. (2010) calculated integral (CIO) and derived the evolution equation for Td as foUowings. 



2 2 



'^d = j^Td - 2jdg{Td - Tg), (Cll) 



where tUpd = \J^T^Q}fnd/md is the dust plasma frequency, and ^dg is the total damping rate due to dust-gas interactions. 
From equation Cll, let denote the excitation coefficient 



= (C12) 



Using typical parameters, we obtain the ratio of sourcc-to-sink from equation (Cll): 



Jdg VO-1/ VlOOe/ VlO^cm 3/ \iych J \ldg J \ rud 

In order to have the increase of dust kinetic energy, the source term must exceed the sink term, i.e., ac/ldg > 1- 

In the low temperature regime, the dust-dust collision is similar to hard-sphere collision. Thus, the variation of Td 
in time is governed by 
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where I = l/{nao) = 1/ {ndTtaj^g) with uhs = 2Aln (Qg/Aer), and A is the screening length. 

Without charge fluctuations, when the collisions between dust particles conserve the energy, the equilibrium temper- 
ature of dust grain is determined by interactions with the ambient gas, so that = Tg. Random charge fluctuations 
provide an additional energy source, and if the coefficient of the first (source) term in the r.h.s. of Eq. (Cll) exceeds 
the damping rate 2r^^, then the dust temperature grows exponentially with time. 

The kinetic coefficient Add is linearly proportional to the maximum scattering angle, which we set equal to unity. 
Hence, we implicitly supposed that there are sufficiently small impact parameters p that ensure scattering at large 
angles, X ^ 1- However, since the lower bound of impact parameters is limited by the particle radius, p ^ a, this 
assumption is only valid if the kinetic temperature is below a certain critical value T^"^. From the relation x ~ Qo/p'^d 
(Lifshitz & Pitaevskii 1981) we readily deduce, 

a e 

so that the actual value of the maximum scattering angle is equal to T" /Td- Therefore, the exponential temperature 
growth described by Eq. (Cll) proceeds until Td reaches the critical value T^'^{^ Ti). At larger temperatures the 
source term in Eq. (Cll) remains constant (which is obtained by replacing Td with T"). The temperature growth 
switches to linear and eventually gets saturated due to the translational damping, which determines the ultimate 
temperature of dust. 

Note that when treating dust-dust collision, wc assume the dust velocity does not change as a result of gas collisions, 
which corresponds to Orbital-Motion-Limited approximation (OML). 

C4^. Analytical estimate for saturated kinetic temperature 

When ac > Jdg, Equation (Cll) indicates that the temperature Td increases exponentially with time. Thus, we call 
this regime instability. When Td becomes larger than the critical temperature T^^ -limit of Coulomb interaction, then 
the source term ac is constant. As Td continues to increase, the sink term increases while the source term remains 
constant. Therefore, the temperature is saturated. From Equation (Cll), we obtain 

2 2 

rrioo _ QpdTjlnrpcr fmr.\ 
Q^Vch ^OL 

Plugging T'^ and parameters z/c/t, cq and Tdn in eq. (C15) it yields 

Ti jin a 0.3(2 -I- z) 

where z ~ Qoe/{akTi) has been used. 
Using the usual expression for thermal velocity 



^ f ik^Td 
V md 

and assuming the MRN dust size distribution, we obtain 

6 X 10^ / T, ^ 



1/2 



(C17) 



where Ui is in cm^'^ and the renormalizing factor a for the damping rate in the subsonic and supersonic regimes is 
given by Eq. (B3) and (B4), respectively. Note that in the latter case a is the function of dust velocity and hence Eq. 
(C18) should be resolved for w^. 

When the source term is smaller than the sink term, ac < Jdg, the dust temperature increases with time until the 
equilibrium with gas is established, and no instability exists. From Equation (Cll) we obtain 



= 1 ^Tn. (C19) 

1 - OLchdn 



As aq increases, TJ"^* increases with oq 
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